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OPTIMIZED SOLUTION OF KEPLER’S EQUATION 


by 

John M. Kohout and Lamar Layton 
Goddard Space Flight Center 


1.0 INTRODUCTION 

1 . 1 General Description of KEPLER 

KEPLER is an IBM 360 computer program used to solve Kepler’s equation for eccentric anomaly: 

E = M + e sin E . 

The double precision input to the program consists of mean anomaly M (in radians) and eccentric- 
ity of the orbit e. The double precision output from program consists of eccentric anomaly E (in 
radians), sin E , and cos E . 

1.2 KEPLER and KEPLR1 

KEPLER has been developed by the authors as a replacement for KEPLR1 , a similar program con- 
tained in the definitive orbit determination system (DODS) used at Goddard Space Flight Center to 
determine orbits for NASA’s scientific satellites. KEPLR1 was not a hastily coded and formulated 
program designed for early replacement as soon as DODS became operational. On the contrary, 
KEPLR1 was developed after a rather extensive research effort that recommended a particular algo- 
rithm, the Myles Standish algorithm.* This algorithm was then implemented by Federal Systems 
Division of IBM under contract to GSFC (see Appendix A). KEPLR1 has been in productive use since 
May 1968. 

KEPLER was developed as a part of an effort to optimize the execution times of frequently 
called subprograms in DODS. The initial thought was simply to recode KEPLR1 in assembly language 
coding (ALC), since, like most other programs in DODS, KEPLR1 was coded in FORTRAN, and 
previous experience with other programs in DODS led the authors to anticipate a 20 to 30 percent 
speed improvement from a FORTRAN-to-ALC recoding. However, a little analysis of the formulation 
used in KEPLR1 led the authors to believe that a new analytic approach to the age-old problem was 


*Cole, Isabella, and Borchers, Raymond V., *‘A Comparison of Some Iterative Techniques for the Solution of Kepler’s Equation”, 
NASA/GSFC Document X-552-67-421, September 1967. 
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in order. As a result, they engaged in a research effort of their own* and developed a completely new 
computer program, KEPLER. 

KEPLER not only solves Kepler’s equation for eccentric anomaly but also outputs accurate values 
for the sine and cosine of the eccentric anomaly. This is important since in almost every case for which 
DODS calls on KEPLR1 to solve Kepler’s equation for E, it then uses E as the input argument to the 
DSIN(X) and DCOS(X) functions. These function calls are unnecessary when KEPLER is used since 
sin E and cos E are part of its output. 

1.3 Outline of Remainder of This Document 

Section 2 of this document is a description of the relative performance of KEPLER versus that of 
KEPLR1. Section 3 is a detailed design and performance description of the newly developed program, 
KEPLER. Section 4 presents the mathematical derivation of the principal formulas used in KEPLER, 
namely the second-order Newton-Raphson differential correction of eccentric anomaly. Section 5 
summarizes the significance of the use of KEPLER and its called module in DODS, DPSC (double 
precision sine/cosine), and recommends the development of related programs. 

Appendix A is a module performance and design description of KEPLR1 prepared for GSFC by 
IBM. Appendixes B and C list the test programs used to compare KEPLER with KEPLR1 . Appendix 
D lists KEPLR1 (in FORTRAN) and KEPLER and its called program, DPSC (in ALC). 

2.0 A COMPARISON OF KEPLER AND KEPLR1 

This section deals with accuracy, speed, core storage requirement, and reentrant properties of the 
KEPLER program. The DODS module KEPLR1 is used as a benchmark for comparison. 

2.1 Accuracy of KEPLER 

Appendix B lists a test program, TEST 1, which is used to exercise both KEPLR1 and KEPLER 
over a full range of the input arguments, M and e. The maximum error produced by each program 
over full ranges of the arguments is obtained by substitution of the solution for eccentric anomaly 
back into Kepler’s equation: 

KEPLR1 error = \E - (M + e sin E)\ = 0.50 X 1(T 15 

KEPLER error = | E - {M + e sin £)| = 0.44 X 10~ 15 . 

2.2 Execution Times 

Appendix C lists a test program, TEST 2, which times the execution of KEPLER and KEPLR1 
on GSFC’s IBM 360/95 computer. Full ranges of e and M are used in this test. The average execution 


*Kohout, J., and Layton, L., “GSFC Optimized Solution of Kepler’s Equation”, NASA/GSFC Document X-541-71-229, May 1971. 
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times reported are based on 400 000 test cases: 

Average KEPLR1 execution time = 314 /is* 
Average KEPLER execution time = 74 /is . 


2.3 Core Requirements 

KEPLR 1 requires 820 bytes of core for itself and 620 bytes of core for its called module, 
IHCLSCN, the standard, release 19 FORTRAN library (FORTLIB) double precision sine/cosine rou- 
tine. KEPLER requires 440 bytes of core for itself and 2360 bytes of core for its called module, 

DPSC. 

2.4 Reentrancy of KEPLER 

Both KEPLER and its called module, DPSC, are reentrant and, therefore, are candidates for the 
high-speed system link pack area. Neither KEPLR1 nor its called module are reentrant. Therefore, 
they may not be stored in the high-speed system link pack area. 

2.5 Optimal Execution Times for KEPLER 

The favorable ratio of execution times reported in Section 2.2 (KEPLR 1 /KEPLER = 314/74 = 
4.24) will be enhanced by factors of 2, 3, and 4, on the IBM 360/95, 360/75, and 360/65, respectively, 
when KEPLER and DPSC are located in the high-speed system link pack area and KEPLR 1 and its 
called module are located in low-speed core. Under these optimal conditions, the execution time ratios 
would be 8.48, 12.72, and 16.96, respectively. 

Furthermore, when KEPLER and DPSC are located in the system link pack area, several concur- 
rent jobs could be calling on KEPLER, and only one copy of KEPLER would be in core. (The core 
storage requirement of KEPLER would be charged to system overhead and not to a particular job.) 

3.0 KEPLER-MODULE PERFORMANCE AND DESIGN DESCRIPTION 

3.1 Language 

KEPLER is written in ALC in order to reduce both execution time and core storage. The core 
storage savings are incidental; the main reason for the use of ALC is to produce a faster executing 
program. It is estimated that the use of ALC for KEPLER is responsible for about one-third of the 
improvement in execution time when that program is used in place of KEPLR 1. 

3.2 Module Size 

KEPLER requires 440 bytes of core storage. 


*This figure includes the time required to calculate sin E and cos E. Without these calculations, the average execution time for 
KEPLR1 is 266 fi s. 
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3.3 Purpose of KEPLER 


KEPLER solves Kepler’s equation for eccentric anomaly when mean anomaly and eccentricity 
are known. The sine and cosine of eccentric anomaly are output by this module. 

3.4 Linkage Information 

3.4.1 Calling Sequence 

KEPLER is invoked via the following call statement: 

CALL KEPLER(MA, ECC, IERR, OUT) , 


where 


MA (input) 
ECC (input) 
IERR (output) 


OUT (output) 
OUT(l) 
OUT(2) 
OUT(3) 


= M, mean anomaly (in radians); 

= e, eccentricity; 

= error code 
= 0, if no error 

= 1, if e is negative or greater than 0.99; 
= 3 word matrix 
= E, eccentric anomaly 
= sin E 
= cos E. 


The qualities M, e, E, sin E, and cos E are double precision floating point numbers. The error code, 
IERR, is a full word integer. 


3.4.2 Called Modules 

KEPLER calls on the reentrant DPSC module. It uses the ALC entry point, SINCOS, which in- 
puts x in FRO and outputs sin x in FRO and cosx in FR2. (Note: FRO is floating point register 0 and 
FR2 is floating point register 2.) 


3.4.3 Calling Modules 

The following modules in DODS call on KEPLER: 

DCCONO, 
CNVRTO, 
UNCALO. 


NEGENO, 

EPTRBO, 

ELCONO, 


3.5 Functional Analysis 

3.5.1 Module Component 1, Main Program 
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3. 5. 1.1 Method 


Kepler’s equation, E = M + e sin E, is solved for E by use of a second-order Newton-Raphson 
iterative algorithm. The derivation of this algorithm is discussed in Section 4 of this document. The 
iterative algorithm is enhanced by four features of KEPLER: 

(1 ) The mean anomaly M is reduced to a value between -it and +ir, and the resultant sign is 
saved. Then, the absolute value of M is used to solve Kepler’s equation. After a solution is obtained, 

E is set equal to 27r - E if the reduced value of M is negative. 

(2) A highly efficient initial estimate algorithm is used to generate E ' , a starting value for the 
iterative process. This algorithm is discussed in Section 3.5.2. 

(3) The SINCOS entry in the DPSC module is used to calculate simultaneously sin E' and cos E' . 

(4) Sum formulas are used to calculate sin ( E ' + C) and cos ( E ' + C) whenever C becomes small 
enough that first- or second-order approximations of sin C and cos C are tolerable. 


3.5.1 .2 Main Program Algorithm 

Step 1 : Error exit if e is negative or greater than 0.99. 

Step 2: Reduce M modulo 2n to range [~7r, +ir] . 

Save sign of M and set M = \M\. 

Step 3: E’ = ESTIMATE {M, e). 

Step 4: Calculate sin E' and cos E' via SINCOS routine. 

Step 5: F = M + e sin E' - E' . 

Step 6: D = 1 - e cos E' . 

Step 7: D' = D + 0.5FesinE'/D. 

Step 8: C= F/D'. 

Step 9: E" = E' + C; store as E' 

Step 10: If |C|> 10~ 5 , return to Step 4. 

Step 11: If |C| < 1CT 8 , skip to Step 1 3. 

Step 1 2: sin E” = sin (E' + C) = (1 - 0.5C 2 ) sin E' + C cos E'; j 

cos E" = cos (£' + C) = (1 - 0.5C 2 ) cos E' - C sin E ' . ) 

Replace sin E' with sin E " , replace cos E' with cos E " , 
and return to Step 5. 

Step 13: sin E" = sin (£' + C) = sin E’ + C cos E'\ I 

cos E" = cos ( E ' + C) - cos E' - C sin E' . ) 

Step 14: If sign of reduced M is positive, skip to Step 16. 

Step 15: Set sin E' = -sin E' and E' = 2ir - E' . 


(See module component 2.) 

(linear correction) 
(first-order derivative) 
(second-order derivative) 
(second-order correction) 
(enhanced E ') 


(second-order sums formulas) 


(first-order sums formulas) 
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Step 1 6: Output E = E' , sin E = sin E", and cos E = cos E", with E, sin E, and cos E accurate to 1 5 

decimal places. 

Step 17: Return to calling program. 

3. 5. 1.3 Explanation of Algorithm 

Step 1 : If e is out of range, no output other than error code is generated. 

Step 2: The reduction of M to the range [ -it, +7r] and the saving of its sign have several advan- 

tages: (1) it increases the precision of the calculation of the linear correction in Step 5 
since M and E' will not exceed it; (2) it simplifies the estimation function (Step 3) since 
M is constrained to the range [ 0, it] ; (3) it provides a convenient test for E being output 
in the range [0, 27 t] (Steps 14-15). 

Step 3: The estimation function is treated in detail in Section 3.5.2. The purpose of this function 

is to provide a sufficiently accurate initial estimate of eccentric anomaly to ensure ( 1 ) 
that the differential correction process defined in Steps 4-9 converges and (2) that this 
convergence takes place in a minimum number of iterations. 

This sequence of steps constitutes one second-order Newton-Raphson iteration. This 
algorithm is considerably more accurate than a first-order differential correction and its 
use has two basic advantages: (1) it converges in fewer iterations than the first-order 
correction and (2) the convergence tolerance e used to terminate the second-order dif- 
ferential correction process may be larger than that for the first-order correction because 
the correction is more accurate. (Other programs use a convergence tolerance of 
5 X 10 -12 , but KEPLER is able to maintain accuracy with a convergence tolerance of 
10~ 8 .) 

If the absolute value of the correction C, is greater than 1 0 -5 , Steps 4-8 are repeated. 

That is, the lengthy SINCOS routine is reexecuted in Step 4 to provide accurate values 
for sin ( E ' + C) and cos ( E ' + C). 

If the absolute value of C is less than 1CT 8 , sufficient convergence is obtained to guaran- 
tee that E is accurate to 1 5 significant digits. Step 1 2 is skipped when this condition is 
met. 

Step 12: When the absolute value of Cis between 10~ 5 and 1 0 — 8 , the algorithm iterates, but the 

lengthy sine/ cosine calculation (Step 4) is replaced by the second-order sum formulas in 
Step 1 2. The largest truncated term in these sum formulas is C 3 /3!. This means that 
when C < 1 0~ 5 , the sum formulas have a relative accuracy of 0. 1 67 X 1 0 -1 5 , which is 
slightly more accurate than the original calculation of sin E' and cos E' in Step 4. 

Step 13: When convergence takes place ( C < 1 0 -8 ), sin E' and cos E' are updated by the first-order 

sum formulas in Step 13. The largest truncated term in the first-order sum formula is 
C 2 / 2!. This means that when C < 1CT 8 , the sum formulas have a relative accuracy of 
0.5 X 10 -16 , which again is more accurate than the original calculation of sin E' and 
cos E' . 


Steps 4-9: 


Step 10: 


Step 1 1 : 
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Step 14: If M is positive after being reduced to the range [- 7 r, + 7 r] , Step 1 5 is skipped 

(sin E > 0, E < 7r). 

Step 15: If 0 > M > -7r, sin £' is set negative, and E f is set equal to 2n - E r . 

Step 16: The quantities E , sin £, and cos £ are sequentially output in a 3 word matrix. 

Step 17: The program is concluded. 

3.5.2 Module Component 2, Initial Estimate 

3.5.2. 1 Method 

As stated in the preceding section, the main purpose of the initial estimate algorithm is to provide 
a starting value for the differential correction process defined in Steps 4-9 of the main program. There 
is an obvious tradeoff between time and accuracy in this initial estimate algorithm. As the initial esti- 
mate is made more accurate, the number of times Steps 4-9 of the main program must be executed is 
reduced. There are, however, constraints on this tradeoff. 

When E 9 is as accurate as one part in 10 5 , the time-consuming sine/cosine calculation step in the 
main program will be executed only once. Therefore, it is desirable to generate E r to this degree of 
accuracy for most combinations of the input parameters e and M. However, it would be uneconomical 
to spend too much time trying to achieve a greater overall accuracy since, regardless of the accuracy 
achieved, the full sine/cosine calculation must be executed at least once in order to further refine E r 
and to generate the sin E and cos E output. 

The algorithm used by KEPLER for the initial estimate was selected only after a large number of 
alternative algorithms were tested and proven to be less efficient.* The name abbreviated Newton- 
Raphson is given to this algorithm because it represents a first-order Newton-Raphson correction in 
truncated precision. It possesses two desirable properties: (1) it is executed in a minimal amount of 
time (37 ps on the IBM 360/75) and (2) for most combinations of M and c, it achieves the desired 
accuracy of one part in 10 5 . 

3.5. 2. 2 Initial Estimate Algorithm (Abbreviated Newton-Raphson) 

The abbreviated Newton-Raphson algorithm consists of two steps. 

Step A: E - M + ez , 0<Af<7r, 

where z is a linear estimate of sin M : 

z = 0.75 M, M < 7r/2; 
z = 0.75 (tt-AO, M> tt/2. 


- , M + e sin E - E 
E ~ E 1 - e cos E 


*Kohou t, J., and Layton, L., “GSFC Optimized Solution of Kepler’s Equation”, NASA/GSFC Document X-541-71-229, May 1971. 
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where sin E and cos E are calculated from the first two terms of the Maclaurin expansions: 


sin x = x - 


3! 


COS X = 1 - — . 

Step B involves several substeps: 

(1 ) Set x = E and S = 1 (S = SWITCH). 

(2) If x > 7t/2, set x = it- x and S = 2. 

(3) If x > 7t/ 4, set x = 7t/2 - x and S = - S. 

(4) Calculate sin E = x - x 2 /6 and cos E = 1 -x 2 12. 

(5) If S is negative, exchange sin E and cos E and set S = -S. 

(6) If S = 2, set cos E = -cosE. 

. - M + e sin E - E 

(7) E' = E+— |n- • 

1 - e cos E 

3. 5. 2. 3 Relation of Component 2 to Main Program 

The initial estimate function (component 2) is linearly coded as Step 3 of the main program 
(component 1 ). Component 2 does not have a separate entry point in KEPLER. 


3.5.3 Flowcharts 

No flowcharts are provided for the main program or the initial estimation function since 
KEPLER’s source code and source code comments directly conform to the logic outlined in Sections 
3. 5. 1.2 and 3.5.2.2. 

3.6 Restrictions and Limitations 

If the input value of e is negative or greater than 0.99, no output other than error code is 
generated. 

3.7 Storage Tables External to KEPLER 
None. 

3.8 Input/Output Device Requirements 
None. 
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3.9 Error Conditions and Recovery 

The error code, IERR, is examined after execution. If IERR = 1 , the input value of e is out of 
range, and hence no other output can be expected. 

3.10 Module Design Technique 

KEPLER is reentrant and, therefore, may be loaded in the system link pack area. KEPLER is 
optimized for fast execution; it is several times faster than existing modules. 

3.1 1 Test Procedures and Results 

The speed and accuracy of KEPLER have been verified by the test programs given in Appen- 
dixes B and C. The results of these tests are summarized in Sections 2.1 and 2.2. 


4.0 SECOND-ORDER NEWTON-RAPHSON METHOD 


Deutsch applies the Newton-Raphson method to the problem of solving Kepler’s equation 
(Reference 1, pp. 24-25). He extends the procedure to include second-order effects, but the final 
equation in the development includes an error in sign, as will be noted later. 


If 


then, 

Let 

Then, 


f(E) = E - M - e sin E , 

/'(£) = 1 -ecosE . 
AE = E, -E 0 . 


AE = 


~(E q - M - e sin E Q ) 
1 - e cos E q 


+ 0[(A£) 2 ] . 


To obtain an expression valid to terms of order (AE) 2 > Deutsch proceeds as follows: 
M = (E q + A E) - e sin ( E Q + A E) 

= E 0 + A E - e(sin E Q cos A E + sin A E cos E Q ) 

(A E-f 


= Eq + A£ - e \ sin E, 


o 


1 - 


+ A E cos E 


o 


then, 


e sin E 


(A E) 2 + ( 1 - e cos E q )A E + (E Q - M - e sin Eq ) - 0 
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Let 


1 


x = 


A E ’ 

A = E q - M - e sin E q , 
B = 1 - e cos Eq , 


C ='- 


e sin E n 


We have then 


Ax 2 + Bx + C = 0 


x = 


_ -B ±yk 2 - 44C 
2A 


Hence, 


AE = 


2A 


-B ±\/b 2 -~4AC 
2A 

-B ± (B - 2AC/B) 


Following Deutsch, we adopt the minus sign as the appropriate choice in the denominator; 
A E= — — 


E q - M - e sin E Q 


-(1 - e cosis 0 ) + (1/2 ){Eq -M-e sin E Q )e sin^/l -e cos E Q ) 


,-l 


AE = 


M - E q + e sin E Q 


1 - e cos Eq + (1/2 ){M - Eq + e sin E 0 )e sin E(1 -e cos E Q ) 


-l 


AE = 


M- Eq + e sin Eq 


1 -efcos£' 0 -(1/2 ){M- Eq + e sini^) sin^/l -e cos E Q ) 1 ] 


(The minus sign within the brackets in the denominator which precedes (1/2) is incorrectly given as a 
plus sign in Reference 1.) 

In general, for functions f(E) for which the relevant derivatives exist, we obtain from a Taylor’s 
series expansion: 

w 

A E= — 

~f'(E 0 ) + f(E Q )f"(E 0 )/2 f\E Q ) ’ 

where terms through (AE) 2 have been included. 
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5.0 CONCLUSIONS AND RECOMMENDATIONS 

5.1 Effect of KEPLER on DODS Performance 

The use of KEPLER in DODS results in a relatively insignificant enhancement of that system 
because, prior to the use of KEPLER, DODS was spending less than three percent of its time solving 
Kepler’s equation. Therefore, even if KEPLER were one hundred times faster than KEPLR1, the 
time saving would not be highly significant in DODS operation. 

5.2 Effect of the DPSC Module on KEPLER 

Of more significance to DODS is the concurrent development of DPSC, with ALC entry points 
SINCOS, DSINX, and DCOSX and FORTRAN entry points, DPSC, DSIN, and DCOS. KEPLER uses 
only the SINCOS entry point. The use of this efficient subroutine to simultaneously calculate sin E' 
and cos E' accounts for about one-third of KEPLER’s enhancement of DODS. Use of ALC and the 
improved mathematical model account for the other two-thirds. 

5.3 Effect of DPSC on DODS 

The inclusion of the DSIN and DCOS entry points in the DPSC module will enhance DODS con- 
siderably more than the inclusion of KEPLER alone. Every sine and cosine calculation in DODS will 
be executed more efficiently since the DSIN and DCOS entry points in DPSC will override the DSIN 
and DCOS entry points in the FORTLIB module, IHCLSCN. 

5.4 Recommended Use of Simultaneous Sine/Cosine Routine 

The DODS formulation contains many situations in which the calculation of both the sine and co- 
sine of a given angle is required. The simultaneous sine/ cosine entry points (DPSC and SINCOS) in the 
DPSC module are now available to DODS programmers who are optimizing DODS modules (such as 
KEPLER) that call for the calculation of both sin (x) and cos (x). The FORTRAN subprogram call, 

CALL DPSC (X, SC) , 

takes about one-half as much time to execute as the separate function calls to IHCLSCN: 

SC(1) = DSIN(X) 

SC(2) = DCOS(X) . 

(When DSIN and DCOS entry points are in the DPSC module, the time enhancement is reduced from 
a factor of 2 to a factor of 1 .5.) 

5.5 Recommended Use of ALC Entry Points in DPSC 

The ALC entry points in the DPSC module enable an ALC program to execute register-to-register 
sine/cosine functions that bypass the highly indirect FORTRAN convention of passing to the function 
program the address of the address of the argument in general register 1 . Besides saving time (the ALC 
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functions are about seven percent faster), the ALC entry points make it possible for some calling pro- 
grams to be written in reentrant code, without using the time-consuming GETMAIN macro. This is 
possible since the ALC functions do not require the storage of an argument list and use only the last 
eight bytes of the save area, which they can share with the calling program. KEPLER is a good exam- 
ple of a second-order reentrant program sharing its save area with the SINCOS routine. 

5.6 Expansion of the DPSC Module 

Because of the frequency of calls to mathematical functions in DODS (and other production pro- 
grams run on GSFC computers), the authors are developing a series of reentrant modules to replace 
the most frequently called function subprogram modules in FORTLIB. The following function sub- 
programs, called TRIGPACK, are all either completed or nearing completion: 


DSQRT(X) , 
DSIN(X) , 
DCOS(X) , 


DTAN(X) , 
DCOT(X) , 
DATAN(X) , 


DATAN2(X,Y) , 
DASIN(X) , 
DACOS(X) . 


The single precision counterparts of these double precision function subprograms are also nearing 
completion. 

Besides the standard FORTRAN entry points, these modules all contain corresponding ALC entry 
points (FORTRAN name with an X suffix— DTANX, for example). The ALC entry points assume that 
the argument is already in floating point register 0. (For the case of the double argument in the 
DATAN2X function, the arguments are assumed to be in floating point registers 0 and 2.) The ALC 
entry points will permit the development of a large number of reentrant second-order subroutines 
since only the last eight bytes of the save area are used and the storage of an argument list is not 
required as a prelude to the subroutines execution. 

5.7 Element Conversion Module 

In conjunction with the development of KEPLER and an optimized reentrant TRIGPACK (Sec- 
tion 5.6), the authors are recoding a DODS module called ELCONO, which contains two inverse sub- 
programs: one for converting position and velocity vectors to osculating Keplerian elements, and the 
other for performing the reverse of this transformation. This third-order module, written in ALC, 
calls on KEPLER and the ALC entry points SINCOS, DSQRTX, DSINX, DATANX, DATAN2X, 
DACOSX, and DTANX in the TRIGPACK modules. It also calls on ALC entry points, VCROSSX, 
VDOTX, and XDOTX in a newly developed vector package. KEPLER and the SINCOS entry in 
DPSC are the heart of this newly optimized module. 
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Appendix A* 


Definitive Orbit Determination System— Model 1 
(Module Performance and Design) 


5. MODULE NAME: KEPLR1— SOLUTION OF KEP LER'S EQUATION 

FOR ECCEN TRIC ANOMALY . 

5.1 LANGUAGE 
FORTRAN IV 

5.2 MODULE SIZE 

The source deck of KEPLR1 consists of 18 executable FORTRAN 
statements and requires 820 bytes of core storage. 

5.3 PURPOSE 

KEPLR1 solves Kepler's equation for eccentric anomaly given mean 
anomaly and eccentricity by the Myles Standish algorithm. 

5.4 INTERFACE INFORMATION 

5.4.1 LINKAGE DEFINITION 

Linkage to this module requires the following CALL statement: CALL 
KEPLR1 (MA, ECC, ERRC, E2) . See Table 1 for the definition of the 
calling sequence arguments. 

5.4.2 INTERFACE BLOCK DIAGRAM 



^Prepared by J. H. Seid, International Business Machines, Inc., under NASA contract NAS5-10022, May 1968. The format employed in 
this appendix is defined in GSFC X-544-70-324: Documentation Standards for the Definitive Orbit Determination System-Per- 
formance and Design Descriptions by R. R. Hohl and Lamar Layton (August 1970). 
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Table 1. Calling Sequence Arguments 


Argument 

Name 

Analytic 

Symbol 

I/O 

Description 

Units 

Format* 

Limits 

Min/Max Dimensions 

MA 

M 

I 

Mean anomaly 

Radians 

LF 

- 

1 

ECC 

e 

I 

Eccentricity 

- 

LF 

0-1 

1 

ERRC 


O 

Error code 

=0, convergence 

<0, no convergence 


LI 


1 

E2 

E 

o 

Eccentric 

anomaly 

Radians 

LF 

- 

1 


♦Format Key 

LF - Long Form Floating Point 
LI - Long Form Integer 


5.4.3 INTERFACE BLOCK DIAGRAM NARRATIVE 
None 

5.4.4 CALLED MODULES 
None 

5.4.5 CALLING MODULES 

ELCONO - Elements Conversion Package 

DC CON 0 - DC Control 
CNVRTO - CONVERT Control 

UNCALO - Unknown Calculation 
CONE GO - CONVERT Normal Equations 
NEGENO - DC Normal Equations 
EPTRBO - EPHEM Tape Record Builder 

5.5 FUNCTIONAL ANALYSIS 

5.5.1 MODULE COMPONENT 1: Solve Kepler's equation 

5. 5. 1.1 Method: 

Given the mean anomaly, M, and the eccentricity, e, the algorithm for 
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computing the eccentric anomaly, E, will be: 

1 . Set error code = 0 

Set limit of number of iterations, MAX = 10 

2 . Set E = 0 

If M =0, go to Step 13 
If M ^0, go to Step 3 

3 . Eq= M + e sin M 

Set number of iterations = 1 

4. F = E Q - (e sin E Q ) - M 

5. D = 1.0 - [ e cos (E q - 0.5F) ] 

6. E = E Q - F/D 

7. If | Eq - E | - TOL <£0, go to Step 13; otherwise continue to Step 8. 

8 . Add 1 to number of iterations 

9. If (number of iterations - MAX) < 0, continue; otherwise go to 
Step 12 

10. E 0 = E 

11 . Return to Step 4 

12 . Set error code = 4 

13. Modulo E by 2 ir 

14 . Return to calling program 

The limit of iterations through Steps 4 to 11 is 10. Thus MAX = 10. If 
this number is exceeded, the error code is set to 4 . 

TOL is the tolerance at which the last significant digit of the difference 
between the previous calculated eccentric anomaly and the present 

_ ^5 

calculated anomaly is allowed. TOL allows an error of Z 5x10 

5. 5. 1.2 Source and Type of Inputs: 

The mean anomaly, M, and the eccentricity, e, will be obtained by this 
module from the calling sequence of the CALL statement which transfers 
control to this function. The format of M and e will be long form 
floating-point . 

5 . 5 . 1 . 3 Destination and Type of Outputs: 

Output of this module will be the eccentric anomaly, E, in long form 
floating-point format. Also the error code, ERRC, in long integer 
format will be returned to the calling sequence of the CALL statement 
which transfers control to this function . 
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Component Level 
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5. 5. 1.5 


Component Level Flowchart Description 

1A - Values for mean anomaly, M, and eccentricity, e, are transferred 
to the module via the calling sequence. The eccentric anomaly, E, 
and the error code, ERRC, are returned via the calling sequence. 
IB - The counter for the number of iterations is set to 1. 

1C - Test to see if mean anomaly is equal to zero. 

ID - Sets the eccentric anomaly, E, equal to zero when the mean 
anomaly, M, is equal to zero and returns to calling module. 

IE - Computes the initial eccentric anomaly, Eq. 

IF - Computes the expression Eq - (e sin Eq) - M using the initial 
computed value of E . 

1G - Computes the expression 1.0 - [e cos (Eq - 0.5F)] using the initial 
computed value of Eq and F computed in IF. 

1H - Computes a more accurate value for eccentric anomaly, 

E = E q - F/D. 

II - Test to see if eccentric anomaly has been determined. If the value 
of the expression |Eq-E | - TOL is equal to or less than zero, E 
has been determined. If E has been determined control is returned 
to the calling module. 

1J - Increase number of iterations by 1 when eccentric anomaly has 
not been determined. 

IK - Test to see if number of iterations is less than or equal to the 
maximum number of iterations allowed. 

1L - If the number of iterations exceeds the maximum, set error code 
equal to 4. Return control to calling module. 

1M - If the number of iterations meets the test, make the initial 

eccentric anomaly, Eq, equal to the computed eccentric anomaly, 

E , and return to 1 F . 

5 . 6 RESTRICTI ONS AND LIMITATIONS 
None 

5 . 7 STORAGE TAB LES EXTERNAL TO MODULE 
None 

5.8 INPUT/OUTPUT DEV ICE R EQ UIREMENTS 
None 


19 



5.9 

5.10 
5.10.1 

5.10.2 

5.10.3 

5.10.4 

5.11 
5.11.1 


ERROR CONDITIONS AND RECOVERY 

Check to see that number of iterations does not exceed the maximum . 

If so, set error code, ERRC =4. 

MODULE DESIGN TECHNIQUE 
MODULARITY REQUIREMENT 
None 

EXPANDABILITY REQUIREMENT 
None 

PARAMETERIZATION 

The maximum number of iterations to determine the eccentric 
anomaly was set equal to 10 . 

SPECIAL FEATURES 
None 

TESTING PROCEDURES AND RESULTS 

UNIT TEST DRIVER DESIGN AND IMPLEMENTATION 

Test 1 - Input data-eccentricity and mean anomaly are read from data 

cards, subroutine KEPLR1 is executed and the results are printed. 

Test 2 - The eccentricity is varied from 0 to 1 by increments of . 05 for 
all values of eccentric anomaly from 0 to 360 degrees, incremented by 
15 degrees in order to compute values of mean anomaly. Subroutine 
KEPLR1 is called for all values of mean anomaly and corresponding 
eccentricity, and the results are printed out. 
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5.11.2 BENCHMARK TESTING 
None 

5.11.3 TEST RESULTS AND ACCURACY 

This item is not covered at the module level; it is included in the system 
evaluation document. 4 '’ 

5.12 GLOSSARY 

A glossary of internal symbols associated with quantities having analytic 
significance is given in Table 2. 

5.13 REFERENCES 

Memorandum from I. Cole to IBM; 1 August 1967 


Table 2. Internal Symbols 


Program 

Symbol 

Analytic 

Symbol 

Description of Term 

Units 

Format 

Limits 
Min/ Max 

El 

E 0 

Eccentricity anomaly 
Compare value 

Radians 

LF 

- 

ITER 

ITER 

Number of iterations 
Completed 

- 

LI 

- 

MAX 

MAX 

Limit of number of 
iterations 

- 

LI 

- 

TOL 

TOL 

Tolerance for 
convergence 

- 

LF 

- 


*Working document in use at GSFC. 
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Appendix B 


Test 1 and Results 


C TEST 1 

C TEST PROGRAM TO TIME TWO KEPLER PROGRAMS 

C PROGRAMMER ANNE BOMFORD 

REAL*8 M,E,EA(3),DE 

500 FORMAT! • 1 ',// ,49X, 'EXECUTION TIMES FOR KEPLER AND KE PLR 1 ' * / / ) 

WRITE! 6*500) 

E=. 10-03 
I T I ME = 2 
CALL I NT I M 0 ( 1 ) 

DO 3 1=1*20 
M= • 1 D— 02 
DO 2 J = 1 * 2 00 
DO 1 K = 1 * 1 00 

1 CALL KEPLER!M,E,IRR,EA> 

2 M = M+ 6 • 28D-02 

3 E = E+ • 49D— 01 

CALL INTIMO! ITIME) 

ITIME=ITIME/ 400000 

506 FORMAT (//* 30X, 'AVERAGE EXECUTION TIME FOR KEPLER ON 360/9 5 ' ♦ IX * I 3 * 

11X* 1 MICROSECONDS' ) 

WRITE!6,506) ITIME 
E=. ID-03 
I T I ME=2 

CALL CrmMO(l) 

DO 6 l=t>?0 
M=. ID-02 
DO 5 J = 1 * 200 
DO 4 K=l, 100 

4 CALL KEPLR1 (M,E* IRR*E2) 

5 M=M*6.28D-02 

6 E = E+.49D-0 1 

CALL INTIMO! ITIME) 

ITIME = ITIM E/400000 

507 FORMAT! //*30X* ‘AVER AGE EXECUTION TIME FOR KE PLR 1 ON 360 /95 ' * 1 X , I 3 * 
llX* 'MICROSECONDS' > 

508 WRITE! 6,507) ITIME 
E = • 1 D-03 

I T I ME =2 

CALL INTIMO! 1 ) 

DO 10 1=1*20 

M= • 1 D-02 

DO 11 J= 1 ,200 

DO 12 K=l, 100 

CALL KEPLRHM,E,IRR*E2) 

E A! 2 ) =DS IN ! E2 ) 

12 EA! 3 ) =DCOS ( E2 ) 

11 M=M+6 • 28D-02 
10 E = E+ • 49D—01 

CALL INTIMO! ITIME ) 

ITIME=ITI ME/ 400000 

607 FORMAT!//, 30X, 'AVERAGE EXECUTION TIME FOR KEPLR1 ON 360/95 ' 1 1 X , I 3 , 

11X,' MICROSECONDS (WITH SINCOS)') 

WR I TE ! 6 , 60 7 ) ITIME 

RETURN 

END 

a*#*:*** i*####*******:##*##**#^ a*## $ * a#t a#c # a** * $ # * ** $ a** * * $ *** * fcjfcajc *Jfe a* a* # 4 = a* * a* 4 : 4 :# 

a* 4=***$ £*4:##*** j{c ^ jJ: j}{ * ^ a*##*##*##**# £ a}: # 4 £ 4 : 4c £ 4c 4c 4e £ * # afc#***#**###*:#* 

EXECUTION TIMES FOR KEPLER AND KEPLRl 
AVERAGE EXECUTION TIME FOR KEPLER ON 360/95 74 MICROSECONDS 

(WITH SINCOS JAVERAGE EXECUTION TIME FOR KEPLRl ON 360/95 314 MICROSECONDS 

(WITHOUT SINCOSJAVERAGE EXECUTION TIME FOR KEPLRl ON 360/95 266 MICROSECONDS 

£ £ £ £#:£ aj=£ 4 ; £ 4= 4= * £ £ £4= 4 : afc ❖ ❖ * 4c:fc # # 4c* * 4= * 4c 4c *4= * ** 4= # 4< 4c*# * £ * * 4c# 4c * £ 4= * 4c4c 4c ** 
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Appendix C 


Test 2 and Results 


C TEST 2 

C TEST PROGRAM TO TEST VALIDITY OF TWO KEPLER PROGRAMS 

C PROGRAMMER ANNE BOMFORD 

R EAL *8 M , E , EA ( 3 ) ,DE 
REAL*8 A , B , ERR » MAX ERR 
REAL*8 E2 

800 FORMAT ( • 1 1 , // , 49X , ‘VALIDITY TEST FOR KEPLER AND KEPLR 1 ' • / / ) 
WRITE<6,8Q0) 

DO 501 K = 1 ,2 

MAXERR=0*0D00 

E = • 1 D-03 

DO 802 1=1,100 

M= . 1 D— 02 

DO 803 J=l,1000 

I F ( K.EQ.2) GO TO 505 

CALL KEPLER(M,E,IRR , E A ) 

ERR=DABS ( EA ( 1 )-M-E*EA< 2) ) 

IFIERR.LE. MAXERR) GO TO 803 
MAXERR = ERR 
A = M 
B = E 

GO TO 803 

505 CALL KEPLR1(M,E,IRR,E2) 

EA ( 1 ) =E2 

EA( 2 ) = DS I N ( E2 ) 

ERR= D AB S ( EA ( 1 )-M-E*EA(2) ) 

IF{ ERR»LE*MAXERR) GO TO 803 
MAXERR =ER R 
A = M 


B=E 

803 M=M+.628D-02 
802 E=E+ • 9800-02 

I F ( K.E0.2) GO TO 808 

806 FORMAT! // ,20X, ‘MAXIMUM ERROR FOR KE P LER= 1 , D 1 0. 3 , • M=‘,E10.3,‘ E = 

1 • , E 1 0 * 3 ) 

WRITE! 6, 806) MAXERR , A, B 
GO TO 501 

807 F0RMAT!//,20X, ‘MAXIMUM ERROR FOR KE PLR 1 = • , 01 0. 3 , • M=»,E10.3,‘ E= 

1 ■ , El 0* 3 ) 

808 WR I T E ! 6 , 8 07 ) MAXERR, A, B 
501 CONTINUE 

RETURN 

END 


4*444444444444444444444:4 44 4444444444444 ^j^#*##:**:*#**:*:*#:*#**##*:*:*#*-^**# 

VALIDITY TEST FOR KEPLER AND KEPLR1 
MAXIMUM ERROR FOR KEPL ER= 0. 4440-15 M= 0.520D 01 E= 0.951D 00 
MAXIMUM ERROR FOR KE P LR1 = 0.500D-15 M= 0.264D 01 E= 0.941D 00 

44 $*******£ s*****#*#****:}:** **4444 4444444 4 44444444444444 4444444 *********** ** ** 
4 4 ************************* ************************************************* 
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Appendix D 


KEPLR1 and KEPLER Source Code 


SUBROUTINE KEPLR1 ( MA, ECC , ERRC , E 2 ) 

THIS SUBROUTINE IS USED TO SOLVE KEPLER'S EQUATION 
FOR ECCENTRIC ANOMALY GIVEN MEAN ANOMALY AND ECCENTRICITY 
BY THE MYLES STANDI SH ALGORITHM. 

IN THE CALLING SEQUENCE MA IS THE MEAN ANOMALY, ECC IS 
THE ECCENTRICITY, ERRC IS THE ERROR CODE FOR NUMBER OF 
ITERATIONS AND ER IS THE ECCENTRIC ANOMALY. 

INTEGER ERRC, MAX, ITER, EC 
DOUBLE PRECISION MA, ECC , El , E2 ,F ,D, ABS 
REAL*8 TOL/. 0 50-10/ , PI 2/6. 28318 5307 179586/ 

MAX= 10 
ERRC=0 
E2=0.0D0 
IF(MA) 3,13,3 
E1=MA+ECC*DSIN(MA) 

I TER= 1 

F = E1— (ECC*DSIN(E1 ) )-MA 
D=1.0D0— ( ECC*DC0S(E1-0.5D0*F ) ) 

E2=E I-F/D 

TEST FOR CONVERGENCE 
I F( DABS (El— E2)— TOL) 13,13,8 
I TER= I TER + 1 

TEST FOR NUMBER OF ITERATIONS 
I F ( ITER-MAX) 10,10,12 
10 E 1 =E2 

GO TO 4 

1 ? PRRf=4 

13 E2= DMOD (E2,PI2) 

I F < E2.LT.0 ) E2 =P 1 2 +E2 

RETURN 

END 


KEPLER CSECT 


* 


ERROR 


USING *, 15 
STM 14,4,12(13) 
LM 1,4, 0(1) 

STEP ONE 

LD 2,0(2) 

LTDR 2,2 

BM ERROR 

CE 2, = E ' . 99 * 

BL *+12 

MV I 3(3), 1 

B EXIT 

MV I 3(3), 0 


SAVE REGISTERS 

ARGUMENT ADDRESSES M» E ,1 ERR , EA 
E 

TEST E FOR MINUS 

TEST E FOR SIZE 
OK IF LESS THAN .99 

* SET ERROR CODE =1 

* AND EXIT 

ZERO ERROR CODE (NO ERROR) 
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* 


* 


QUAD1 


0CT1 


STEP TWO 


LD 

4,0(1) 


M HOLD FOR SIGN 

LPDR 

0 ,4 


M ABSOLUTE 

CD 

0, TWOP I 


OK IF LESS THAN 2PI 

BL 

*+20 


* 

LDR 

6,0 


* OTHERWISE 

DD 

6 ,TWOPI 


* REDUCE M 

AW 

6,ZER014 


* TO MODULO 2PI 

MD 

6, TWOP I 


* 

SDR 

0,6 


* 

CD 

0 ,PI 


* 

BL 

*+12 


* IF X GREATER THAN PI 

LCDR 

4,4 


* AND REVERSE SIGN 

LCDR 

0,0 


* AND SET X = 2PI-X 

AD 

0 ,TWOPI 


* 

STE 

4,52 ( 13) 


SIGN OF M SAVED 

STD 0,56(13) 

EP THREE-A 


M ABSOLUTE BETWEEN 0 AND PI 

CE 

0 ,*E' 1.57* 

* 


BL 

*+10 

* 


LCDR 

0,0 

* 

QUICK AND DIRTY SIN(M) 

AD 

0,PI 

* 


HDR 

0,0 

* 

= ( 3/4 ) X X= M M BETWEEN 

HDR 

4,0 

* 

OR X = PI — M M BETWEEN 

ADR 

0 ,4 

* 


MER 

0,2 


E*SIN(E ) 

AD 

0 ,56( 13) 

EA 

= M + E*SIN (M ) , 1ST ESTIMATE 

STD 

0,0(4) 

HOLD AS EA (POSITIVE) 


O-PI/2 

PI/2-PI 


STEP THREE-B 


LA 

CE 

BL 

LA 

LCDR 

AO 

CE 

BL 

LCR 

LCDR 

AE 

LCDR 

HOR 

MER 

MER 

ME 

AD 

ADR 

LTR 

BP 

LCR 

LDR 

LDR 

LDR 

BCT 

LCDR 


S= 1 (ASSUMES EA IN QUAD 1 ) 


0,1 

0 ,=E* 1.5707963* 

QUAD1 OK 1ST QUADRANT 

0,2 OTHERWISE SET 

0,0 AND 

0 ,PI EA= P I -EA 

0,=E' .7853981* 

0CT1 OK 1ST 

0,0 
0,0 
0, =E 1 
2,0 
4,0 
2,4 
4,2 

4, =E' .6666667* 

2, ONEX 


S=2 


EA LESS THAN P I /2 


OCTANT 

* OTHERWISE 

* AND 


SET S - * 


1.5707963* * 

-X 


EA = PI/2 -EA 


* 


EA LESS PI/4 


0 ,4 
0,0 
* + 12 
0,0 

4.2 
2,0 
0 ,4 
0 , *+6 

2.2 


TEST 

OK 1ST 
* 

* 

* 

* 

* —COS ( EA ) 

* -COS(EA) 


X/2 
-X2/2 
-X3/4 
-X3/6 
1-X2/2 
X-X3/6 
S FOR + 

OCTANT 
OTHERWISE 
AND 

EXCHANGE 
SIN(EA) AND 
- IN 


= COS(EA) 
= SIN(EA) 
OR - 


SET S =+l OR +2 


+ IN 


COS(EA) 
QUAD 1 
QUAD 2 
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ME 

0,012) 

E*SIN( EA ) 


ME 

2,0(2) 

— E*COS(EA) 


AD 

2, ONEX 

1— E*COS( EA ) =0 1ST ORDER 


AD 

0 ,56( 13) 

M+E*SIN(EA> 


SD 

0,0(4) 

M+E*SIN(EA> -EA =F 


DER 

0 ,2 

F/D =C 1ST ORDER 


AD 

0,0(4) 

EA+C 


STD 

0,0(4) 

STORE AS EA ESTIMATE (POSITIVE) 

♦ 

STEP FOUR 



L 

15,=V( SINCOS) 


BALR 

14*15 



USING 

*, 14 


♦ 

STEP FIVE 


ITERATE STD 

0 ,8(4) 

SIN(EA) 


MD 

0,0(2) 

E*SIN(EA) 


HDR 

4,0 

.5*E*SIN(EA) HOLD 


AD 

0,56 ( 13) 

M+E*SIN( EA) 


SD 

0 ,0(4) 

M+E*SIN(EA> -EA =F 

* 

STEP SIX 




LCDR 

6,2 

-COS(EA) 


HD 

6,0(2) 

— E*COS (E A ) 


AD 

6, ONEX 

1-E*C0S(EA) =D 1ST ORDER 

* 

STEP SEVEN 



MER 

4,0 

F*.5*E*SIN(EA) 


DER 

4,6 

F*.5*E*SIN(EA)/D 


ADR 

6,4 

D +F*.5*E*SIN( EA)/D = D 2ND ORDER 

♦ 

STEP EIGHT 



DDR 

0 *6 

F/D = C 2ND ORDER 


LDR 

4,0 

SAVE C 

* 

STEP NINE 



AD 

0,0(4) 

EA+C 


STD 

0 ,0(4) 

SAVE AS ENHANCED EA 

* 

STEP TEN 




LPER 

6,4 

C ABSOLUTE 


CE 

6»=E • l.E-5 * * RETURN FOR FULL SINCOS AND ITERATE 


BCR 

2, 15 

* IF C GREATER THAN .00001 

* 

STEP ELEVEN 



CE 

6,=E'1.E 

-8' * CONVERGENCE 


BL 

OUT 

* WHERE C LESS THAN .00000001 

* 

STEP TWELVE 



LCDR 

0,4 

-C 2ND ORDER SUMS FORMULAE 


HDR 

6,4 

C/2 C BETWEEN 10**-5, 10**-8 


MDR 

0,6 

-C*C/2 


AD 

0 ,ONEX 

I— C*C/2 = COS(C) 


LDR 

6,4 

C = SIN (C ) 


MDR 

6,2 

SIN ( C )*COS ( EA) 


MDR 

2,0 

COS (C ) *COS( EA ) 


MD 

4,8(4) 

SIN ( C )*SIN ( EA) 


MD 

0,8(4) 

COS (C )* SI N( EA ) 


ADR 

0 ,6 

SIN( EA+C )*SI N( C )*COS( EA) +COS (C)*SIN(EA> 


SDR 

2,4 

COS ( EA+C ) =COS ( C ) *COS ( EA )-S I N( C ) * SI N( EA ) 


BR 

14 

ITERATE 

* 

STEP THIRTEEN 


OUT 

LDR 

0 ,4 

C 1ST ORDER SUMS FORMULAE * 


MDR 

0,2 

C*COS!EA) C LESS THAN 10**-8 * 


* 

* 

* 

* 

* 

* 

* 

♦ 

* 

* 

* 
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AD 

0 ,8(4) 

SIN( EA+C ) 

= SIN(EA) +C*COS ( EA 


MD 

4,8(4) 

C*SIN ( EA ) 



SDR 

2,4 

COS( EA+C) 

= COS (EA) -C*S IN ( EA 

* 

STEP FOURTEEN 




TM 

52(13), 

12 8 

CHECK SIGN OF M 


BZ 

PLUS 


OK IF M POSITIVE 

* 

STEP FIFTEEN 




LCDR 

0,0 

OTHERWISE 

COMPLIMENT SIN(EA) 


LO 

4 »TWOPI 

* AND 



SD 

4,0(4) 

* SET 



STD 

4,0(4) 

* EA = 2PI-EA 

* 

STEP SIXTEEN 



PLUS 

STD 

0 ,8(4) 

OUTPUT SIN(EA) 


STD 

2,16(4) 

OUTPUT COS( EA ) 


* STEP SEVENTEEN 


EXIT 

LM 

14,4,12(13) RESTORE 

REGISTERS 


BR 

14 AND RETURN 

TWOPI 

DC 

D* 6. 283 18 5307 179 58 64 * 

2PI 

PI 

DC 

X '413243F6A8885A31* 

PI 

ZER014 

DC 

X • 4E 00000000000000 ' 


ONEX 

DC 

END 

X '40FFFFFFFFFFFFFF • 

=E*.99* 

=E ' 1 .57 • 

=E'l. 5707963* 

=E. 7853981 • 

=E'. 6666667* 

=V ( S I NCOS) 

= E ' 1 .E-5 * 



* 

* 

* 


I NT I MO 


=E 1 1 • E— 8 1 
START 0 

BC 15,12(15) 

DC X ■ 7 1 
DC CL7 1 INTIMO 1 
STM 14,12,12(13) 
BALR 12,0 
USING * , 12 
LR 3,1 


BRANCH AROUND CONSTANTS 
ESTABLISH A HALF-WORD BOUNDARY 
NAME 

SAVE THE REGISTERS 
BASE REGISTER 


INTERV 


L 7,0(3) 

CLI 3 ( 7 ) , X • 2 1 

BC 8, INTERV 

STIMER T ASK , TUI NT VL= INTER 

LM 14,12,12(13) 

MV I 12 ( 1 3 ) , X 1 FF • 

BCR 15,14 

L 5, INTER 

TTIMER 


CHECK THE INDICATOR 
BRANCH IF SECOND ENTRY 
SET THE TIMER 
RESTORE THE RGGISTERS 
INDICATE CONTROL R ETURNEC 
RETURN TO CALLING PROGRAM 
LOAD MAXIMUM TIME 


SR 

M 

ST 
LM 
MV I 
BCR 
DS 

INTER DC 

TWSIX DC 

END 


5,0 

4, TWSIX 
5,0(7) 

14,12, 12( 13) 
12(13) , X 1 F F • 
15,14 
OF 

F 1 2147483647 • 
F # 26* 


DETERMINE ELAPSED TIME 
CONVERT UNITS TO MICRO SECONDS 
STORE TIME IN RETURN LOCATION 
RESTORE REGI STERS 
INDICATE CONRROL RETURNED 
RETURN TO CALLING PROGRAM 
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DPSC 


SI NCOS 


CSECT 




USING 

*, 15 



STM 

14,2, 12(13) 

SAVE REGS 

LM 

I » 2 » 0 ( 1 ) 


X,SC ADDRESSES 

LD 

0»0< I) 


X TO FRO 

LA 

1 5 » S I NCOS 


SINCOS ADD 

BALR 

14,15 


SIN TO FRO, COS TO FRO 

STD 

0 ,0(2) 


OUTPUT SIN(X) 

STD 

2,8(2) 


OUTPUT COS(X) 

LM 

14,2, 12( 13) 

RESTORE REGS 

BR 

14 


AND RETURN 

ENTRY 

S INCOS 



USING 

*,15 



LDR 

6,0 


SIGN OF SIN(A) 

LPDR 

0,0 


X ABSOLUTE 

CD 

0 t TWOPI 

❖ 


BL 

*+20 

* 

IF X GREATER THAN 2PI 

LDR 

2,0 

* 

REDUCE 

DD 

2t TWOP I 

* 

TO 

AW 

2 , ZERO 14 

* 

MODULO 

MD 

2, TWOP I 

* 

TWOPI 

SDR 

0 ,2 

* 


CD 

0,PI 


* 

BL 

* + 12 


* IF X GREATER THAN PI 

LCDR 

6,6 


* REVERSE SIGN OF SINE 

LCDR 

0 ,0 


* AND SET X = 2PI-X 

AD 

0, TWOPI 


* 

LA 

0 ,2 

SET SWITCH = 2 (COS +, NO SII 

CD 

0, P I 0VER2 


* 

BL 

*+14 


* IF X GREATER THAN PI /2 

LA 

0,1 


* SET X = P I-X 

LCDR 

0 ,0 


* AND SET SWITCH =1 (COS -) 

AD 

0,PI 


* 

CD 

0 , P I0VER4 


* IF X GREATER THAN P I /4 

BL 

*+16 


* SET X = PI/4 +PI /4 -X 

LCR 

0,0 


* AND SET SWITCH - (S 

LCDR 

0,0 


* 

AD 

0 , P 1 0VER4 



AD 

0, P I 0VER4 


* 

LER 

2,0 



AU 

2 , =X ’45000000 ' 4500000X 1ST HEX DIGIT ( 

STE 

2 , 68 ( 13) 



L 

1 ,68 ( 13) 


45000001 

SLL 

1 .3 


81 = CONSTANT INDEX 1 = 

LDR 

2,0 


X TO FR 2 AND FRO 

MDR 

2*, 2 


X*X =X2 

LDR 

4,2 


X2 

MD 

4,S5( 1) 



AD 

4, S4 ( 1) 



MDR 

4,2 



AD 

4,S3( 1) 



MDR 

4,2 



AD 

4, S2 ( 1) 



MDR 

4,2 



AD 

4, SI ( 1) 



MDR 

0 ,4 


SIN(X) VIA 9TH DEGREE OD 


= INDEX 


1 = 0 , 1 , 2,... ,12 


LDR 

4,2 

X2 

CE 

4 , =E • • 25 ' 

X2 VS 1/4 

BL 

*+18 

8TH DEGREE COS IF X LESS THAN 1/2 

MD 

2 »C6— 64 ( 1) 


AD 

2»C5 ( 1) 

10TH DEGREE COS IF X GREATER THAN 1/2 

MDR 

2,4 


B 

*+8 


MD 

2,C5( 1) 


AD 

2,C4( 1) 


MDR 

2,4 


AD 

2»C3( 1) 


MDR 

2,4 


AD 

2,C2 ( 1) 


MDR 

2,4 


AD 

2,C1 ( 1) 

COStXT VIA 8 OR 10 DEGREE POLYNOMIAL 

LTR 

0 ,0 

TEST SWITCH FOR SIGN 

BP 

*+12 

NO EXCHANGE IF + 

LCR 

0,0 

OTHERWISE SET SWITCH +1 OR +2 

LDR 

4,2 

AND 

LDR 

2,0 

EXCHANGE 

LDR 

0,4 

SIN AND COS 

BCT 

0 ,*+6 

* IF SWITCH = 1 

LCDR 

2,2 

* SET COSINE - 

LTDR 

6,6 

* IF SIGN OF SIN + 

BCR 

10,14 

* EXIT 

LCDR 

0,0 

OTHERWISE SET SIGN - 

BR 

14 

AND EXIT 


USING *, 15 
ENTRY DCOS,DCOSX 
DCOS L 1,0(1) 

LD 0,0(1) 

LNDR 0,0 

AD 0»P I0VER4 

AD 0 »PI0VER4 

LA 15, DS INX 
BR 15 

USING *,15 
DCOSX LNDR 0,0 

AD 0,P I0VER4 

AD 0 »PI0VER4 

LA 15, DS INX 

BR 15 

ENTRY DSIN,DS INX 
OSIN L 1,0(1) 

LD 0,0(1) 

LA 15,12(15) 

USING *,15 

DS 1 NX LDR 6,0 

LPDR 0,0 

CD 0 »TWOP I 

BL *+20 

LDR 2,0 

DD 2, TWOP I 

AW 2 »ZER014 

MD 2,TW0PI 
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coss 


SDR 0,2 

CD 0,PI 

BL *+10 

LCDR 6,6 

SD 0,PI 

CD 0»P I0VER2 

BL *+10 

LCDR 0,0 

AD 0,PI 

CD 0»P I0VER4 

BNL COSS 

LER 2,0 

AU 2 ,=X '45000000* 

STE 2,68(13) 

L 1,68(13) 

SLL 1,3 

LDR 2,0 

MDR 2,2 

LDR 4,2 

MD 4,S5( 1) 

AD 4 » S4 ( 1 ) 

MDR 4,2 

AD 4 , S3 ( 1 ) 

MDR 4,2 

AD 4 , S2 ( 1 ) 

MDR 4,2 

AD 4 , SI ( 1 ) 

MDR 0,4 

LTDR 6,6 

BCR 10,14 

LCDR 0,0 

BR 14 

LCDR 0,0 

AD 0»P I0VER4 

AD 0 , P I 0VER4 

LER 2,0 

AU 2 , =X * 45000000 • 

STE 2,68(13) 

L 1,68(13) 

SLL 1,3 

MDR 0,0 

LDR 2,0 

CE 2 , =E * . 25 ' 

BL *+18 

MD 0 ,C6— 64 ( 1 ) 

AD 0,C5( 1) 

MDR 0,2 

B *+8 

MD 0 ,C5 ( 1 ) 

AD 0,C4( 1 ) 

MDR 0,2 

AD 0,C 3 ( 1) 

MDR 0,2 

AD 0,C2 ( 1 ) 

MDR 0,2 
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AD 

0,CI< 1) 





LTDR 

6 » 6 





BCR 

10t 14 





LCDR 

0 ,0 





BR 

14 





DS 

OD 




ZER014 

DC 

X*4E00000000000000' 




TWOPI 

DC 

0 '6.2831853071795864* 

2PI 



PI 

DC 

X'413243F6A8885A31 ' 

PI 



PI0VER2 

DC 

X '411921FB54442D18' 

PI/2 


PI0VER4 

DC 

X'40C90FDAA22168C2 ' 




* 

SINE 

COEFFICIENTS 




SI 

DC 

X' 40FFFFFFFFFFFFFA • 


0.0 



DC 

X '40FFFFFFFFFFFF37* 


0.6250000000D- 

-01 


DC 

X ' 40FFFFFFFFFFFD7C * 


0.12 50000 000 D 

00 


DC 

X ' 40FFFFFFFFFFE8D8 ' 


0.18 750000 OOD 

00 


DC 

X ' 40FFFFFFFFFFC96A ' 


0 • 2500000000D 

00 


DC 

X '40FFFFFFFFFF49BA' 


0. 3125000000D 

00 


DC 

X'40FFFFFFFFFCFB23' 


0*37 50000 OOOD 

00 


DC 

X '40FFFFFFFFF01C57 • 


0.437 5000000D 

00 


DC 

X'40FFFFFFFFCC5C90» 


0 • 5000000000D 

00 


DC 

X '40FFFFFFFF608C98 • 


0. 5625000000D 

00 


DC 

X'40FFFFFFFE63DA15' 


0.6250000 OOOD 

00 


DC 

X '40FFFFFFFC1E8C81* 


0 .68 75000000D 

00 


DC 

X'40FFFFFFF85A9960 ' 


0.7 5 00000 OOOD 

00 

S2 

DC 

X • C02AAAAAAAAA7 13F • 


0.0 



DC 

X'C02AAAAAAAA931F5 ' 


0.62 50000 000 D- 

-01 


DC 

X ' C02 AA AAAAAA8FD27 * 


0.12500000 OOD 

00 


DC 

X'C02AAAAAAAA2 D215 ' 


0. 1875000000D 

00 


DC 

X 'C02 AAAAAAA9E6858 * 


0. 2500000000D 

00 


DC 

X • C02AAAAAAA8C DOO A • 


0*31 25000 OOOD 

00 


DC 

X ' C02AAAAAAA4CBBA7' 


0.37500000 OOD 

00 


DC 

X'C02AAAAAA947FA52 ' 


0*4375000 OOOD 

00 


DC 

X • C02AAAAAA7 IE 40EB ' 


0.50000000 OOD 

00 


DC 

X'C02AAAAAA2826B13' 


0.5625 000 OOOD 

00 


DC 

X ' C02 AAAAA9802504 8 ' 


0. 625 000 00 OOD 

00 


DC 

X' C02AAAAA84FD178E • 


0*6875 000 OOOD 

00 


DC 

X 'C02AAAAA69F0419B* 


0. 7 5000 000 OOD 

00 

S3 

DC 

X'3F222221FB15B5CF • 


0.0 



DC 

X '3F22222212475F04* 


0.6250000000D- 

■01 


DC 

X'3F2222221B88AD9A' 


0.1250000000D 

00 


DC 

X '3F2222221225BEF1' 


0.18750000 OOD 

00 


DC 

X'3F2222221106E925' 


0. 2 5 0000 000 OD 

00 


DC 

X '3F22222 203229637* 


0.31250000 OOD 

00 


DC 

X' 3F222221D93FF490 • 


0.37 50000 OOOD 

00 


DC 

X ' 3F22222 15CBB1C5D* 


0.4375000000D 

00 


DC 

X'3F 22222093 4E 194 1 * 


0 • 5 000 000 OOOD 

00 


DC 

X • 3F2 222 1F3A705D7F • 


0.56250000 OOD 

00 


DC 

X * 3F22221CBB9904E F * 


0.6250000000D 

00 


DC 

X'3F222218FE37316B* 


0.68750000 OOD 

00 


DC 

X'3F2222146F331054* 


0.7 500000000 D 

00 

S4 

DC 

X ' BDCFBE6DAB3 3F92D* 


0.0 



DC 

X* B0D008956DAAE A4B ' 


0*625000000 OD- 

*01 


DC 

X 'BDD00C4AC2C8225A' 


0.12500000 OOD 

00 


DC 

X'BDD00C1538OE D130' 


0 . 18 75000000D 

00 
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o # 


DC X 1 BDD00C4960797744 1 

DC X* BDD00C008 1D853F7 ' 

DC X • BDD00B3D87649F44' 

DC X 1 BDD0099428EEBB 20 1 

DC X • BDD00789E59C8C04' 

DC X 1 BDD004BB8 BF 343 6D 

DC X • BDD0008 16BFF2313 1 

DC X 1 BDCFF B44 C6 16A71 4 

DC X *BDC C F5CCED99F6BD' 

S5 DC X*BD27ABC0019B0975 

DC X 1 3C2 7402 0A16 A8780 1 

DC X * 3C2DC 3E3 A8 A96C9 1 

DC X 1 3C2 DE2C5F202B5DE 1 

DC X * 3C2E01B434E997D5 

DC X * 3C2DF9F7766C5869 

DC X * 3C2DE4A5BC08 1B7 1 

DC X •3C2DC260BE20B1AB 1 

DC X * 3C2DA291 A650688 7 

DC X * 3C2D7F6C 3453E649 

DC X • 3C2D545C68C 1566F 

DC X 9 3C2D284AE4CA5DCB’ 

DC X* 3C2D00E6F7D092D6 

COSINE COEFFICIENTS 
I DC X* 41 10000000000000 

DC X *4110000000000000 

DC X* 41 10000000000010 

DC X * 40FFFFFFFFFFF6B5 

DC X 1 40FFFFFF FFFF02CC 

DC X * 40FFFFFFFFF9D8B0 

DC X* 40FFFFFFFFDB4BF3 

DC X * 40FFFFFFFF66CA3C 

DC X* 40FFFFFFFFBFB694 

DC X *411 000000002779F 

DC X • 41 1 000000000C48 4 

DC X *40FFFFFFFFF0ACD4 

DC X • 40FFF FFF FF328278 

C2 DC X *C0800000000008FC 

DC X * C08 000000000 1 F56 

DC X * C08 000000 00 0A56E 

DC X * C07FFFFF FFFB4 22C 

DC X *C07FFFFFFFC18410 

DC X*C07FFFFFFEF2E9A8 

DC X * C 07 FFF FFF BA 00 34 A 

DC X * C07FF FFF F2 4D4282 

DC X *C07FFFFFFB7C7615 

DC X * CO 80000002 24 OB 38 

DC X • C0800000006D 49C3 

DC X* C07FFFFFFF 1 A5797 

DC X *C 07FFFF FF88 B5 A64 

C3 DC X* 3FAA AAAAACFAD91B 

DC X • 3F AAAAAAABB4F60A 

DC X * 3FA AAAAAAC8F2FE5 

DC X * 3F AAAAAA9BD824DF 

DC X 1 3FAA AAAA487F9F90 

DC X * 3FAAA AA985B06FCA 


0.2500000000D 00 
0*31250000000 00 
0*37500000000 00 
0*43750000000 00 
0.5000000000D 00 
0*5625 000 00 OD 00 
0*62500000000 00 
0*68750000000 00 
0* 75000 OOOOOD 00 
0*0 

0.6250000000D-01 
0 • 1 2 50000 000 D 00 
0. 18 75000000D 00 
0* 2500000000D 00 
0. 3125000000D 00 
0 *37 50000 00 OD 00 
0*437 5000000D 00 
0*50000000000 00 
0*562 50 OOOOOD 00 
0*62 50000 OOOD 00 
0*68 75000000D 00 
0*7 500000000 D 00 

0.0 

0* 62 5 00 OOOOOD- 01 
0* 1250000000D 00 
0*18 750000000 00 
0* 2 500000 OOOD 00 
0*31250000000 00 
0* 37 50000 OOOD 00 
0 • 437 5000000D 00 
0 *50000000000 00 
0*56250000000 00 
0*62500000000 00 
0 • 68 75000000D 00 
0*7 5 00000 OOOD 00 
0.0 

0 *62 500000000—0 1 
0* 1250000000D 00 
0. 1875000000D 00 
0* 250 00 OOOOOD 00 
0*3125 000 OOOD 00 
0.3750000000D 00 
0*4375000 OOOD 00 
0* 5000000000D 00 
0*5625 000 00 OD 00 
0* 6250000000D 00 
0 * 68 75 000 OOOD 00 
0* 750000 0000 D 00 
0.0 

0. 62 5000 00 OOD— 01 
0 * 12 50000 OOOD 00 
0 • 18 75000000D 00 
0* 2 5 00000 OOOD 00 
0* 31 2 50 OOOOOD 00 
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DC 

X* 3FAAAAA756C2 D236 • 

0. 37 50000000 D 

010 

DC 

X ’ 3FAAAAA2D85E052B ' 

0*43750000000 

00 

DC 

X*3FAAAAA8A34E1B3E 1 

0.5000000000D 

00 

DC 

X • 3FAAAAAB6650E57D' 

0. 56250000 OOD 

00 

DC 

X • 3FAAAAAABBBF861 4' 

0.6250000000D 

00 

DC 

X • 3FAAAAAA5544824F • 

0.68 75000000D 

00 

DC 

X* 3FAAAAA8E1C6081B ' 

0.7 500000000D 

00 

DC 

X • BE5B05E544ACFAED* 

0.0 


DC 

X*BE5B05B3A3D232D3* 

0*62500000000 

-01 

DC 

X •BE5B05B13C670BD8* 

0. 1250000000D 

00 

DC 

X • BE5B059A0ECC D046 * 

0.1875000000D 

00 

DC 

X ' BE5B 0563A45C8 E91 • 

0.25000000 OOD 

00 

DC 

X*BE5B0511816D0F42 ' 

0.3125000000D 

00 

DC 

X •BE5B046D420B0F0B' 

0. 3750000000D 

00 

DC 

X* BE5B0375D64C08AC * 

0.4375000 00 OD 

00 

DC 

X 'BE5B053B6D3 50E39* 

0. 50000000 OOD 

00 

DC 

X*BE5B05CF9CE1EFA5* 

0. 5625000000D 

00 

DC 

X ‘BE5B05AF4705E015* 

0.62500000 OOD 

00 

DC 

X , BE5B059FF606305B ' 

0.6875000000D 

00 

DC 

X * BE5B0576C215745A* 

0.75000000 OOD 

00 

DC 

X • 301B83EA0F58760F * 

0.0 


DC 

X *301 A0381E55BED90* 

0.6250000000D- 

-01 

DC 

X* 3D19FD902FF6A01A * 

0. 1 2 50000000 D 

00 

DC 

X • 3D19F150857FB484* 

0.18 750000 OOD 

00 

DC 

X*3D19E3CB6B016FA1* 

0. 2500000000D 

00 

DC 

X '3D19D6B96D105CAD* 

0. 3125000000D 

00 

DC 

X*3D19C48042D 34407* 

0.37 50000000D 

00 

DC 

X '3D19B076C277FE83* 

0.437 5000000D 

00 

DC 

X*3D19F45FAA3F5386 • 

0. 5000000000D 

00 

DC 

X » 3D1 A040A698 23422 • 

0.562 5000000D 

00 

DC 

X* 3D1 A0107B2 4 39E 18 ' 

0 .62500000000 

00 

DC 

X '3D19FFE547902218* 

0.68 750000 OOD 

00 

DC 

X* 3D19FD9A C88121EB* 

0.7 500000000D 

00 

DC 

X * BB4010A69645256D* 

0. 5000000000D 

00 

DC 

X'BB4A BF02911 31076 * 

0.5625000000D 

00 

OC 

X • BB48F9D4FC6AC452 * 

0. 6250000000D 

00 

DC 

X*BB487121A5139E61* 

0.6875000 00 OD 

00 

DC 

X • BB479FD89AA019C2 * 

0. 750000 OOOOD 

00 

END 

=X *45000000* 




=E • • 2 5 * 
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